using Distributions
using StatsPlots
using LaTeXStrings
using CSV
using DataFrames
using StatisticalRethinking
using Logging
using Turing
using MCMCChains

# setting default attributes for plots
default(label=false)

4.1 Why normal distributions are normal

Code 4.1

n = rand(Uniform(-1, 1), 1000, 16);
pos = sum.(eachrow(n));
density(pos)

Code 4.2

prod(1 .+ rand(Uniform(0, .1), 12))

2.0638427726255317

Code 4.3

u = Uniform(0, .1)
growth = prod.(eachrow(1 .+ rand(u, 10000, 12)));

density(growth; label="density")
# overlay normal distribution
μ = mean(growth)
σ = std(growth)
plot!(Normal(μ, σ); label="Normal")

Code 4.4

big = prod.(eachrow(1 .+ rand(Uniform(0, 0.5), 10000, 12)));
small = prod.(eachrow(1 .+ rand(Uniform(0, 0.01), 10000, 12)));
density([big, small]; layout=(2, 1))

Code 4.5

density(log.(big))

4.2 A language for describing models

Code 4.6

w = 6
n = 9
p_grid = range(0, 1; length=100)
bin_dens = [pdf(Binomial(n, p), w) for p in p_grid]
uni_dens = [pdf(Uniform(0, 1), p) for p in p_grid];
posterior = bin_dens .* uni_dens
posterior /= sum(posterior);

4.3 Gaussian model of height

Code 4.7

d = DataFrame(CSV.File("data/Howell1.csv"));

Code 4.8

describe(d)

4 rows × 7 columns

	variable	mean	min	median	max	nmissing	eltype
	Symbol	Float64	Real	Float64	Real	Int64	DataType
1	height	138.264	53.975	148.59	179.07	0	Float64
2	weight	35.6106	4.25242	40.0578	62.9926	0	Float64
3	age	29.3444	0.0	27.0	88.0	0	Float64
4	male	0.472426	0	0.0	1	0	Int64

Code 4.9

precis(d)

┌────────┬────────────────────────────────────────────────────────────┐
│  param │    mean      std     5.5%      50%    94.5%      histogram │
├────────┼────────────────────────────────────────────────────────────┤
│ height │ 138.264  27.6024  81.1086   148.59  165.735  ▁▁▁▂▂▂▂▂▂██▆▁ │
│ weight │ 35.6106  14.7192   9.3607  40.0578  54.5029  ▁▃▄▄▃▂▃▆██▅▃▁ │
│    age │ 29.3444  20.7469      1.0     27.0   66.135      █▆▆▆▆▃▃▁▁ │
│   male │  0.4724   0.4997      0.0      0.0      1.0    █▁▁▁▁▁▁▁▁▁█ │
└────────┴────────────────────────────────────────────────────────────┘

Code 4.10

d.height

544-element Vector{Float64}:
 151.765
 139.7
 136.525
 156.845
 145.415
 163.83
 149.225
 168.91
 147.955
 165.1
 154.305
 151.13
 144.78
   ⋮
 156.21
 152.4
 162.56
 114.935
  67.945
 142.875
  76.835
 145.415
 162.56
 156.21
  71.12
 158.75

Code 4.11

d2 = d[d.age .>= 18,:];

Code 4.12

plot(Normal(178, 20); xlim=(100, 250))

Code 4.13

plot(Uniform(0, 50), xlim=(-10, 60), ylim=(0, 0.03))

Code 4.14

size = 10_000
sample_μ = rand(Normal(178, 20), size)
sample_σ = rand(Uniform(0, 50), size);
prior_h = [rand(Normal(μ, σ)) for (μ, σ) in zip(sample_μ, sample_σ)];

p1 = density(sample_μ; title="μ")
p2 = density(sample_σ; title="σ")
p3 = density(prior_h; title="prior_h")

plot(p1, p2, p3, layout=@layout [a b; c])

Code 4.15

sample_μ = rand(Normal(178, 100), size)
prior_h = [rand(Normal(μ, σ)) for (μ, σ) in zip(sample_μ, sample_σ)];

density(prior_h)
vline!([0, 272])

Code 4.16

μ_list = range(150, 160; length=100)
σ_list = range(7, 9; length=100)

log_likelihood = [
    sum(logpdf(Normal(μ, σ), d2.height))
    for μ ∈ μ_list, σ ∈ σ_list
]
log_prod = log_likelihood .+ [
    logpdf(Normal(178, 20), μ) + logpdf(Uniform(0, 50), σ)
    for μ ∈ μ_list, σ ∈ σ_list
];

max_prod = maximum(log_prod)
prob = @. exp(log_prod - max_prod);

Code 4.17

# note the transposition, that's due to Julia matrix order
contour(μ_list, σ_list, prob')

Code 4.18

heatmap(μ_list, σ_list, prob')

Code 4.19

indices = collect(Iterators.product(1:length(μ_list), 1:length(σ_list)));
sample_idx = wsample(vec(indices), vec(prob), 10_000; replace=true)
sample_μ = μ_list[first.(sample_idx)]
sample_σ = σ_list[last.(sample_idx)];

Code 4.20

scatter(sample_μ, sample_σ; alpha=0.1)

Code 4.21

p1 = density(sample_μ)
p2 = density(sample_σ)
plot(p1, p2, layout=(2,1))

Code 4.22

println(hpdi(sample_μ, alpha=0.11))
println(hpdi(sample_σ, alpha=0.11))

[153.83838383838383, 155.15151515151516]
[7.262626262626263, 8.191919191919192]

Code 4.23

d3 = sample(d2.height, 20);

Code 4.24

μ_list = range(150, 170; length=200)
σ_list = range(4, 20; length=200)

log_likelihood = [
    sum(logpdf(Normal(μ, σ), d3))
    for μ ∈ μ_list, σ ∈ σ_list
]
log_prod = log_likelihood .+ [
    logpdf(Normal(178, 20), μ) + logpdf(Uniform(0, 50), σ)
    for μ ∈ μ_list, σ ∈ σ_list
]

max_prod = maximum(log_prod)
prob2 = @. exp(log_prod - max_prod)

indices = collect(Iterators.product(1:length(μ_list), 1:length(σ_list)));
sample2_idx = wsample(vec(indices), vec(prob2), 10_000; replace=true)
sample2_μ = μ_list[first.(sample2_idx)]
sample2_σ = σ_list[last.(sample2_idx)]

scatter(sample2_μ, sample2_σ; alpha=0.1)

Code 4.25

density(sample2_σ)
μ = mean(sample2_σ)
σ = std(sample2_σ)
plot!(Normal(μ, σ); label="Normal")

Code 4.26

d = DataFrame(CSV.File("data/Howell1.csv"));
d2 = d[d.age .>= 18,:];

Code 4.27

@model function model_height(height)
    μ ~ Normal(178, 20)
    σ ~ Uniform(0, 50)
    height ~ Normal(μ, σ)
end

model_height (generic function with 1 method)

Code 4.28

Logging.disable_logging(Logging.Warn)
m4_1 = sample(model_height(d2.height), NUTS(), 1000)
print(MCMCChains.header(m4_1))

Iterations        = 501:1:1500
Number of chains  = 1
Samples per chain = 1000
Wall duration     = 7.03 seconds
Compute duration  = 7.03 seconds
parameters        = μ, σ
internals         = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size

Code 4.29

display.(describe(m4_1; q=[0.055, 0.945]));

Summary Statistics
  parameters       mean       std   naive_se      mcse        ess      rhat    ⋯
      Symbol    Float64   Float64    Float64   Float64    Float64   Float64    ⋯

           μ   154.6030    0.4265     0.0135    0.0191   870.3501    1.0003    ⋯
           σ     7.7681    0.2812     0.0089    0.0107   824.1419    0.9998    ⋯
                                                                1 column omitted

Quantiles
  parameters       5.5%      94.5% 
      Symbol    Float64    Float64 

           μ   153.9036   155.2826
           σ     7.3407     8.2574

Code 4.30

init_vals = [mean(d2.height), std(d2.height)]
chain = sample(model_height(d2.height), NUTS(), 1000, init_theta=init_vals)

Chains MCMC chain (1000×14×1 Array{Float64, 3}):

Iterations        = 501:1:1500
Number of chains  = 1
Samples per chain = 1000
Wall duration     = 0.67 seconds
Compute duration  = 0.67 seconds
parameters        = μ, σ
internals         = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size

Summary Statistics
  parameters       mean       std   naive_se      mcse         ess      rhat   ⋯
      Symbol    Float64   Float64    Float64   Float64     Float64   Float64   ⋯

           μ   154.6031    0.4326     0.0137    0.0109   1043.6194    0.9991   ⋯
           σ     7.7640    0.3119     0.0099    0.0069   1043.3823    1.0007   ⋯
                                                                1 column omitted

Quantiles
  parameters       2.5%      25.0%      50.0%      75.0%      97.5% 
      Symbol    Float64    Float64    Float64    Float64    Float64 

           μ   153.7697   154.3036   154.6012   154.8870   155.5128
           σ     7.2058     7.5439     7.7487     7.9656     8.4026

Code 4.31

@model function model_height(height)
    μ ~ Normal(178, 0.1)
    σ ~ Uniform(0, 50)
    height ~ Normal(μ, σ)
end

m4_2 = sample(model_height(d2.height), NUTS(), 1000)
display.(describe(m4_2; q=[0.055, 0.945]));

Summary Statistics
  parameters       mean       std   naive_se      mcse         ess      rhat   ⋯
      Symbol    Float64   Float64    Float64   Float64     Float64   Float64   ⋯

           μ   177.8616    0.1005     0.0032    0.0023   1009.0981    1.0009   ⋯
           σ    24.6179    0.9806     0.0310    0.0193    813.0156    0.9990   ⋯
                                                                1 column omitted

Quantiles
  parameters       5.5%      94.5% 
      Symbol    Float64    Float64 

           μ   177.7031   178.0197
           σ    23.0864    26.2440

Code 4.32

cov(hcat(m4_1[:μ], m4_1[:σ]))

2×2 Matrix{Float64}:
  0.181935   -0.0033216
 -0.0033216   0.0790992

Code 4.33

c = cov(hcat(m4_1[:μ], m4_1[:σ]))
cov2cor(c, diag(c))

2×2 Matrix{Float64}:
  1.0       -0.230813
 -0.230813   1.0

Code 4.34

# resetrange is needed due to bug in MCMCChains: https://github.com/TuringLang/MCMCChains.jl/issues/324
# once it will be fixed, direct sampling from the chain will be possible
samp_chain = sample(resetrange(m4_1), 10_000)
samp_df = DataFrame(samp_chain)
first(samp_df, 5)

5 rows × 2 columns

	μ	σ
	Float64	Float64
1	154.69	7.25339
2	154.35	8.26152
3	154.326	8.1032
4	155.242	7.61533
5	154.405	7.23033

Code 4.35

precis(samp_df)

┌───────┬─────────────────────────────────────────────────────────┐
│ param │    mean     std     5.5%      50%    94.5%    histogram │
├───────┼─────────────────────────────────────────────────────────┤
│     μ │ 154.602  0.4221  153.911  154.607  155.266       ▁▂▆█▄▁ │
│     σ │  7.7721  0.2798   7.3453   7.7561   8.2567  ▁▁▂▆█▇▄▂▁▁▁ │
└───────┴─────────────────────────────────────────────────────────┘

Code 4.36

data = hcat(m4_1[:μ], m4_1[:σ])
μ = mean(data, dims=1)
σ = cov(data)
mvn = MvNormal(vec(μ), σ)
post = rand(mvn, 10_000);
print(mvn)

FullNormal(
dim: 2
μ: [154.60295741552895, 7.768086720261678]
Σ: [0.18193494174975705 -0.0033216008639284323; -0.0033216008639284323 0.07909916601423655]
)

4.4 Linear predictions

Code 4.37

d = DataFrame(CSV.File("data/Howell1.csv"));
d2 = d[d.age .>= 18,:];
# fancy way of doing scatter(d2.weight, d2.height)
@df d2 scatter(:weight, :height)

Code 4.38

Random.seed!(2971)
N = 100
a = rand(Normal(178, 20), N)
b = rand(Normal(0, 10), N);

Code 4.39

p = hline([0, 272]; ylims=(-100, 400), xlabel="weight", ylabel="hegiht")
title!(L"\beta \sim \mathcal{N}(\mu=0,\sigma=10)")

x_mean = mean(d2.weight)
xlims = extrema(d2.weight)  # getting min and max in one pass

for (α, β) ∈ zip(a, b)
    plot!(x -> α + β * (x - x_mean); xlims=xlims, c=:black, alpha=0.3)
end
display(p)

Code 4.40

b = rand(LogNormal(0, 1), 10_000)
density(b, xlims=(0, 5), bandwidth=0.1)

Code 4.41

Random.seed!(2971)
N = 100
a = rand(Normal(178, 20), N)
b = rand(LogNormal(0, 1), N);

Code 4.42

d = DataFrame(CSV.File("data/Howell1.csv"));
d2 = d[d.age .>= 18,:]
xbar = mean(d2.weight)

@model function height_regr_model(weight, height)
    a ~ Normal(178, 20)
    b ~ LogNormal(0, 1)
    μ = @. a + b * (weight - xbar)
    σ ~ Uniform(0, 50)
    height ~ MvNormal(μ, σ)
end

m4_3 = sample(height_regr_model(d2.weight, d2.height), NUTS(), 1000)
m4_3 = resetrange(m4_3);

Code 4.43

@model function height_regr_model_exp(weight, height)
    a ~ Normal(178, 20)
    log_b ~ Normal(0, 1)
    μ = @. a + exp(log_b) * (weight - xbar)
    σ ~ Uniform(0, 50)
    height ~ MvNormal(μ, σ)
end

m4_3b = sample(height_regr_model_exp(d2.weight, d2.height), NUTS(), 1000);

Code 4.44

m4_3_df = DataFrame(m4_3)
precis(m4_3_df)

┌───────┬────────────────────────────────────────────────────────┐
│ param │    mean     std     5.5%      50%    94.5%   histogram │
├───────┼────────────────────────────────────────────────────────┤
│     a │ 154.606  0.2625  154.175  154.616  155.014  ▁▁▂▅██▅▂▁▁ │
│     b │  0.9014  0.0423    0.835   0.9021   0.9665      ▁▃██▃▁ │
│     σ │  5.1071  0.1978   4.8019   5.0921   5.4553     ▁▁▅█▅▂▁ │
└───────┴────────────────────────────────────────────────────────┘

Code 4.45

round.(cov(Matrix(m4_3_df)), digits=3)

3×3 Matrix{Float64}:
  0.069  -0.0    -0.001
 -0.0     0.002  -0.0
 -0.001  -0.0     0.039

Code 4.46

p = @df d2 scatter(:weight, :height; alpha=0.3)

chain = resetrange(m4_3)
samples = sample(chain, 1000)

a_map = mean(samples[:a])
b_map = mean(samples[:b])
plot!(x -> a_map + b_map*(x-xbar))

Code 4.47

post = sample(m4_3, 1000)
post_df = DataFrame(post)
post_df[1:5,:]

5 rows × 3 columns

	a	b	σ
	Float64	Float64	Float64
1	155.241	0.939497	5.32578
2	154.435	0.890103	5.15068
3	154.393	0.95766	5.12624
4	154.644	0.975131	5.02818
5	154.719	0.833828	5.12795

Code 4.48

N = 10
dN = d2[1:N,:]

@model function height_regr_model_N(weight, height)
    a ~ Normal(178, 20)
    b ~ LogNormal(0, 1)
    m_weight = mean(weight)
    μ = @. a + b * (weight - m_weight)
    σ ~ Uniform(0, 50)
    height ~ MvNormal(μ, σ)
end

mN = sample(height_regr_model(dN.weight, dN.height), NUTS(), 1000)
mN = resetrange(mN);

Code 4.49

post = sample(mN, 20)
post_df = DataFrame(post);

xlims = extrema(d2.weight)
ylims = extrema(d2.height)
p = @df dN scatter(:weight, :height; xlims=xlims, ylims=ylims)
title!("N = $N"; xlab="weight", ylab="height")

x_mean = mean(dN.weight)
for (a, b) ∈ zip(post_df.a, post_df.b)
    plot!(x -> a + b * (x-x_mean); c="black", alpha=0.3)
end
display(p)

Code 4.50

post = sample(m4_3, 1000)
post_df = DataFrame(post)
μ_at_50 = @. post_df.a + post_df.b * (50 - xbar);

Code 4.51

density(μ_at_50; lw=2, xlab=L"\mu|weight=50")

Code 4.52

PI(μ_at_50; prob=0.89)

2-element Vector{Float64}:
 158.5792439935818
 159.65866731906047

Code 4.53

μ = StatisticalRethinking.link(post_df, [:a :b], d2.weight, xbar);
μ = hcat(μ...);
Base.size(μ), μ[1:5,1]

((1000, 352), [157.56293534189655, 157.70353697789378, 156.97988672219105, 156.85965609569527, 157.2632797224988])

Code 4.54

weight_seq = 25:70
μ = StatisticalRethinking.link(post_df, [:a :b], weight_seq, xbar);
μ = hcat(μ...);
Base.size(μ), μ[1:5,1]

((1000, 46), [137.04072934028633, 136.64103095564798, 135.80904017330352, 136.73820023564727, 137.81265964505803])

Code 4.55

p = plot()
for i in 1:100
    scatter!(weight_seq, μ[i,:]; c=:blue, alpha=0.2)
end
display(p)

Code 4.56

μ_mean = mean.(eachcol(μ))
μ_PI = PI.(eachcol(μ))
μ_PI = vcat(μ_PI'...);

Code 4.57

@df d2 scatter(:weight, :height; alpha=0.2, xlab="weight", ylab="height")
plot!(weight_seq, [μ_mean μ_mean]; c=:black, fillrange=μ_PI, fillalpha=0.3)

Code 4.58

post = sample(m4_3, 1000)
post = DataFrame(post)

weight_seq = 25:70
μ = map(w -> post.a + post.b * (w - xbar), weight_seq)
μ = hcat(μ...)
μ_mean = mean.(eachcol(μ))
μ_CI = PI.(eachcol(μ));

Code 4.59

sim_height = simulate(post, [:a, :b, :σ], weight_seq .- xbar);
Base.size(sim_height), sim_height[1:5,1]

((1000, 46), [131.0600330744535, 131.07613382456364, 144.0054107126664, 133.0809884065339, 139.1993923582497])

Code 4.60

height_PI = PI.(eachcol(sim_height))
height_PI = vcat(height_PI'...);

Code 4.61

@df d2 scatter(:weight, :height; alpha=0.2, xlab="weight", ylab="height")
plot!(weight_seq, [μ_mean μ_mean]; c=:black, fillrange=μ_PI, fillalpha=0.3)
plot!(weight_seq, [μ_mean μ_mean]; c=:black, fillrange=height_PI, fillalpha=0.3)

Code 4.62

post = sample(m4_3, 10_000)
post = DataFrame(post)
sim_height = simulate(post, [:a, :b, :σ], weight_seq .- xbar)
height_PI = PI.(eachcol(sim_height))
height_PI = vcat(height_PI'...);

@df d2 scatter(:weight, :height; alpha=0.2, xlab="weight", ylab="height")
plot!(weight_seq, [μ_mean μ_mean]; c=:black, fillrange=μ_PI, fillalpha=0.3)
plot!(weight_seq, [μ_mean μ_mean]; c=:black, fillrange=height_PI, fillalpha=0.3)

Code 4.63

post = sample(m4_3, 1000)
post = DataFrame(post)

sim_height = [
    [
        rand(Normal(a + b * (w - xbar), σ))
        for (a, b, σ) ∈ zip(post.a, post.b, post.σ)
    ]
    for w ∈ weight_seq
]
sim_height = hcat(sim_height...)

height_PI = PI.(eachcol(sim_height));
height_PI = vcat(height_PI'...);

@df d2 scatter(:weight, :height; alpha=0.2, xlab="weight", ylab="height")
plot!(weight_seq, [μ_mean μ_mean]; c=:black, fillrange=μ_PI, fillalpha=0.3)
plot!(weight_seq, [μ_mean μ_mean]; c=:black, fillrange=height_PI, fillalpha=0.3)

4.5 Curves from lines

Code 4.64

d = DataFrame(CSV.File("data/Howell1.csv"))
scatter(d.weight, d.height; alpha=0.3)

Code 4.65

d[!, :weight_s] = standardize(ZScoreTransform, d.weight)
d[!, :weight_s2] = d.weight_s.^2;

@model function height_regr_m2(weight_s, weight_s2, height)
    a ~ Normal(178, 20)
    b1 ~ LogNormal(0, 1)
    b2 ~ Normal(0, 1)
    μ = @. a + b1 * weight_s + b2 * weight_s2
    σ ~ Uniform(0, 50)
    height ~ MvNormal(μ, σ)
end

m4_5 = sample(height_regr_m2(d.weight_s, d.weight_s2, d.height), NUTS(), 1000)
m4_5 = resetrange(m4_5)
m4_5_df = DataFrame(m4_5);

Code 4.66

precis(m4_5_df)

┌───────┬───────────────────────────────────────────────────────────┐
│ param │    mean     std     5.5%      50%    94.5%      histogram │
├───────┼───────────────────────────────────────────────────────────┤
│     a │ 146.054  0.3497  145.501  146.059  146.576  ▁▁▂▃▅██▇▄▂▁▁▁ │
│    b1 │ 21.7228  0.2941  21.2498  21.7259  22.1891    ▁▁▁▄▇█▇▄▂▁▁ │
│    b2 │ -7.8033  0.2589  -8.1943  -7.8236  -7.3754      ▁▁▂▄█▆▅▂▁ │
│     σ │  5.8137   0.179   5.5276   5.8055   6.1068   ▁▁▁▃▆██▆▄▂▁▁ │
└───────┴───────────────────────────────────────────────────────────┘

Code 4.67

Random.seed!(1)
df = sample(m4_5_df, 1000)
weight_seq = range(-2.2, 2; length=30)

# explicit logic of link
mu = [
    df.a + df.b1 * w_s + df.b2 * w_s^2
    for w_s ∈ weight_seq
]

mu = hcat(mu...)
mu_mean = mean.(eachcol(mu))
mu_PI = PI.(eachcol(mu))
mu_PI = vcat(mu_PI'...)

# explicit logic of sim
sim_height = [
    [
        rand(Normal(row.a + row.b1 * w_s + row.b2 * w_s^2, row.σ))
        for row ∈ eachrow(df)
    ]
    for w_s ∈ weight_seq
]
sim_height = hcat(sim_height...);

height_PI = PI.(eachcol(sim_height))
height_PI = vcat(height_PI'...);

Code 4.68

p_square = @df d scatter(:weight_s, :height; alpha=0.3, title="Square poly")
plot!(weight_seq, mu_mean; c=:black)
plot!(weight_seq, [mu_mean mu_mean]; c=:black, fillrange=mu_PI, fillalpha=0.3)
plot!(weight_seq, [mu_mean mu_mean]; c=:black, fillrange=height_PI, fillalpha=0.3)

Code 4.69

d[!, :weight_s3] = d.weight_s.^3;

@model function height_regr_m3(weight_s, weight_s2, weight_s3, height)
    a ~ Normal(178, 20)
    b1 ~ LogNormal(0, 1)
    b2 ~ Normal(0, 1)
    b3 ~ Normal(0, 1)
    μ = @. a + b1 * weight_s + b2 * weight_s2 + b3 * weight_s3
    σ ~ Uniform(0, 50)
    height ~ MvNormal(μ, σ)
end

m4_6 = sample(height_regr_m3(d.weight_s, d.weight_s2, d.weight_s3, d.height), NUTS(), 1000)
m4_6 = resetrange(m4_6)
m4_6_df = DataFrame(m4_6)
precis(m4_6_df)

┌───────┬───────────────────────────────────────────────────────────┐
│ param │    mean     std     5.5%      50%    94.5%      histogram │
├───────┼───────────────────────────────────────────────────────────┤
│     a │ 146.399  0.3158  145.888  146.405  146.898    ▁▁▁▃▅██▆▃▁▁ │
│    b1 │ 15.2346  0.4765  14.4594  15.2243  15.9857       ▁▁▁▄█▅▁▁ │
│    b2 │  -6.213  0.2561  -6.6352  -6.2127  -5.8004     ▁▁▂▅██▄▂▁▁ │
│    b3 │  3.5764  0.2239   3.2084    3.584   3.9364       ▁▂▄██▃▁▁ │
│     σ │  4.8605   0.147   4.6456   4.8546   5.0868  ▁▁▁▄▇█▇▄▂▁▁▁▁ │
└───────┴───────────────────────────────────────────────────────────┘

Random.seed!(1)
df = sample(m4_6_df, 1000)
weight_seq = range(-2.2, 2; length=30)

# explicit logic of link
mu = [
    df.a + df.b1 * w_s + df.b2 * w_s^2 + df.b3 * w_s^3
    for w_s ∈ weight_seq
]

mu = hcat(mu...)
mu_mean = mean.(eachcol(mu))
mu_PI = PI.(eachcol(mu))
mu_PI = vcat(mu_PI'...)

# explicit logic of sim
sim_height = [
    [
        rand(Normal(row.a + row.b1 * w_s + row.b2 * w_s^2 + row.b3 * w_s^3, row.σ))
        for row ∈ eachrow(df)
    ]
    for w_s ∈ weight_seq
]
sim_height = hcat(sim_height...);

height_PI = PI.(eachcol(sim_height))
height_PI = vcat(height_PI'...);

p_cubic = @df d scatter(:weight_s, :height; alpha=0.3, title="Cubic poly")
plot!(weight_seq, mu_mean; c=:black)
plot!(weight_seq, [mu_mean mu_mean]; c=:black, fillrange=mu_PI, fillalpha=0.3)
plot!(weight_seq, [mu_mean mu_mean]; c=:black, fillrange=height_PI, fillalpha=0.3)

plot(p_square, p_cubic; layout=(1, 2))

Code 4.70 and 4.71

Looks like Julia plots don't support change of ticks proposed in the book. But much more natural way will be to remap values we're plotting back to the original scale. Example of this is below.

μ = mean(d.weight)
σ = std(d.weight)
w = @. d.weight_s * σ + μ
scatter(w, d.height; alpha=0.3)

w_s = @. weight_seq * σ + μ
plot!(w_s, mu_mean; c=:black)
plot!(w_s, [mu_mean mu_mean]; c=:black, fillrange=mu_PI, fillalpha=0.3)
plot!(w_s, [mu_mean mu_mean]; c=:black, fillrange=height_PI, fillalpha=0.3)

Code 4.72

d = DataFrame(CSV.File("data/cherry_blossoms.csv", missingstring="NA"))
precis(d)

┌────────────┬─────────────────────────────────────────────────────────┐
│      param │   mean      std    5.5%     50%    94.5%      histogram │
├────────────┼─────────────────────────────────────────────────────────┤
│       year │ 1408.0  350.885  867.77  1408.0  1948.23  ████████████▂ │
│        doy │ 104.54    6.407   94.43   105.0    115.0       ▁▂▆██▅▂▁ │
│       temp │ 6.1419   0.6636    5.15     6.1   7.2947       ▁▄▆█▄▂▁▁ │
│ temp_upper │ 7.1852   0.9929  5.8977    7.04   8.9023       ▂██▃▁▁▁▁ │
│ temp_lower │ 5.0989   0.8503  3.7876   5.145     6.37       ▁▁▁▂▇█▂▁ │
└────────────┴─────────────────────────────────────────────────────────┘

Code 4.73

d2 = d[completecases(d[!,[:doy]]),:]
d2 = disallowmissing(d2[!,[:year,:doy]])
num_knots = 15
knots_list = quantile(d2.year, range(0, 1; length=num_knots));

Code 4.74

using BSplines

basis = BSplineBasis(3, knots_list)

16-element BSplineBasis{Vector{Float64}}:
 order: 3
 breakpoints: [812.0, 1036.0, 1174.0, 1269.0, 1377.0, 1454.0, 1518.0, 1583.0, 1650.0, 1714.0, 1774.0, 1833.0, 1893.0, 1956.0, 2015.0]

Code 4.75

p1 = plot(basis)
scatter!(knots_list, repeat([1], num_knots); xlab="year", ylab="basis", legend=false)

Code 4.76

This way of calucalting bsplines is slightly slower, than shown in the book (with pre-calculated matrix) but it is much cleaner in my perspective.

You can do comparison yourself by precalculating spline matrix outside of the model and do matrix multiplication in the model instead of spline evialutaion. Example of doing this is at code block 4.79

@model function model_splines(year, doy)
    w ~ MvNormal(zeros(length(basis)), 1)
    a ~ Normal(100, 10)
    s = Spline(basis, w)
    μ = a .+ s.(year)
    σ ~ Exponential(1)
    doy ~ MvNormal(μ, σ)
end

m4_7 = sample(model_splines(d2.year, d2.doy), NUTS(0.65; init_ϵ = 9.765625e-5), 1000)

Chains MCMC chain (1000×30×1 Array{Float64, 3}):

Iterations        = 501:1:1500
Number of chains  = 1
Samples per chain = 1000
Wall duration     = 13.95 seconds
Compute duration  = 13.95 seconds
parameters        = w[1], w[13], σ, w[3], w[12], w[8], w[10], w[14], w[11], w[9], w[16], a, w[2], w[4], w[6], w[5], w[7], w[15]
internals         = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size

Summary Statistics
  parameters       mean       std   naive_se      mcse         ess      rhat   ⋯
      Symbol    Float64   Float64    Float64   Float64     Float64   Float64   ⋯

        w[1]    -0.6812    0.9172     0.0290    0.0129   1813.8789    0.9991   ⋯
        w[2]    -0.8555    0.8626     0.0273    0.0229   2705.8346    0.9990   ⋯
        w[3]     0.8252    0.8027     0.0254    0.0263   1626.8872    1.0002   ⋯
        w[4]     0.6955    0.7619     0.0241    0.0170   1808.3373    0.9990   ⋯
        w[5]     0.5630    0.7818     0.0247    0.0151   2370.8672    1.0008   ⋯
        w[6]    -1.2257    0.7915     0.0250    0.0182   2674.2308    0.9994   ⋯
        w[7]     0.2788    0.7483     0.0237    0.0179   2438.4965    0.9996   ⋯
        w[8]     1.8091    0.7712     0.0244    0.0159   2157.9901    1.0000   ⋯
        w[9]    -0.0526    0.7338     0.0232    0.0140   2400.1936    0.9991   ⋯
       w[10]     1.9439    0.8019     0.0254    0.0167   2163.3095    0.9995   ⋯
       w[11]     0.6797    0.7354     0.0233    0.0217   1598.7700    1.0000   ⋯
       w[12]     1.2764    0.7565     0.0239    0.0163   2692.3364    0.9990   ⋯
       w[13]     1.2634    0.7793     0.0246    0.0155   2460.8566    0.9994   ⋯
       w[14]    -0.7788    0.8162     0.0258    0.0147   2129.6796    0.9998   ⋯
       w[15]    -2.9017    0.8942     0.0283    0.0160   2867.7978    0.9991   ⋯
       w[16]    -2.7443    0.8905     0.0282    0.0183   2674.3391    0.9994   ⋯
           a   104.2685    0.3292     0.0104    0.0131   1152.1117    1.0032   ⋯
           σ     6.1107    0.1658     0.0052    0.0036   2781.4362    0.9991   ⋯
                                                                1 column omitted

Quantiles
  parameters       2.5%      25.0%      50.0%      75.0%      97.5% 
      Symbol    Float64    Float64    Float64    Float64    Float64 

        w[1]    -2.5031    -1.2788    -0.6836    -0.0650     1.0794
        w[2]    -2.5142    -1.4874    -0.8544    -0.2671     0.8429
        w[3]    -0.7560     0.2729     0.8494     1.3813     2.3222
        w[4]    -0.7929     0.1645     0.6935     1.1964     2.1724
        w[5]    -0.9255     0.0087     0.5661     1.0865     2.1346
        w[6]    -2.8200    -1.7612    -1.2422    -0.6835     0.3608
        w[7]    -1.2176    -0.2135     0.2979     0.7884     1.7322
        w[8]     0.2887     1.2753     1.8034     2.3243     3.3166
        w[9]    -1.4817    -0.5690    -0.0725     0.4282     1.3626
       w[10]     0.3610     1.4107     1.9837     2.5043     3.4097
       w[11]    -0.8514     0.1983     0.6660     1.1726     2.1750
       w[12]    -0.1970     0.7800     1.2584     1.8213     2.7546
       w[13]    -0.3652     0.7719     1.2357     1.8412     2.7361
       w[14]    -2.3775    -1.3282    -0.7610    -0.2397     0.8342
       w[15]    -4.7354    -3.4840    -2.8894    -2.3133    -1.1709
       w[16]    -4.5349    -3.3456    -2.7131    -2.1420    -1.0208
           a   103.6152   104.0409   104.2619   104.4908   104.8871
           σ     5.7975     5.9989     6.1047     6.2182     6.4409

Code 4.77

post = DataFrame(m4_7)

# convert columns w[*] into single column w
w_df = DataFrames.select(post, r"w")
post = DataFrames.select(post, Not(r"w"))
post[!,:w] = Vector.(eachrow(w_df))

# vector of 16 average w values
w_mean = mean.(eachcol(w_df))
p2 = plot(basis .* w_mean)
scatter!(knots_list, repeat([1], num_knots); xlab="year", ylab="basis × weight")

Code 4.78

# explicit link logic
μ = [
    row.a .+ Spline(basis, row.w).(d2.year)
    for row ∈ eachrow(post)
]
μ = hcat(μ...);

μ_PI = PI.(eachrow(μ))
μ_PI = vcat(μ_PI'...);

p3 = @df d2 scatter(:year, :doy; alpha=0.3)
μ_mean = mean.(eachrow(μ_PI))
plot!(d2.year, [μ_mean, μ_mean]; c=:black, fillrange=μ_PI, fillalpha=0.3, alpha=0)

plot(p1, p2, p3; layout=(3, 1))

Code 4.79

How to build the model with explicit spline matrix calculation

basis = BSplineBasis(3, knots_list)

# list of splines with 1 only at corresponding basis index
splines = [
    Spline(basis, [float(idx == knot) for idx ∈ 1:length(basis)])
    for knot ∈ 1:length(basis)
]

# calculate each spline for every year. Resulting matrix B is 827x16
B = [
    map(s -> s(year), splines)
    for year in d2.year
]
B = vcat(B'...);


# do not need years parameter anymore, all the information is in B matrix
@model function model_splines_matrix(doy)
    w ~ MvNormal(zeros(length(basis)), 1)
    a ~ Normal(100, 10)
    μ = a .+ B * w
    σ ~ Exponential(1)
    doy ~ MvNormal(μ, σ)
end

m4_7alt = sample(model_splines_matrix(d2.doy), NUTS(0.65; init_ϵ = 0.0001953125), 1000)

Chains MCMC chain (1000×30×1 Array{Float64, 3}):

Iterations        = 501:1:1500
Number of chains  = 1
Samples per chain = 1000
Wall duration     = 4.91 seconds
Compute duration  = 4.91 seconds
parameters        = w[1], w[13], σ, w[3], w[12], w[8], w[10], w[14], w[11], w[9], w[16], a, w[2], w[4], w[6], w[5], w[7], w[15]
internals         = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size

Summary Statistics
  parameters       mean       std   naive_se      mcse         ess      rhat   ⋯
      Symbol    Float64   Float64    Float64   Float64     Float64   Float64   ⋯

        w[1]    -0.6977    0.9024     0.0285    0.0118   2200.5259    0.9990   ⋯
        w[2]    -0.8587    0.8423     0.0266    0.0204   2112.9445    1.0009   ⋯
        w[3]     0.8351    0.7801     0.0247    0.0201   2085.0565    0.9990   ⋯
        w[4]     0.6791    0.6900     0.0218    0.0234   1563.5078    0.9990   ⋯
        w[5]     0.5912    0.8079     0.0255    0.0150   2520.6431    0.9990   ⋯
        w[6]    -1.2373    0.7507     0.0237    0.0177   2291.9338    0.9995   ⋯
        w[7]     0.2550    0.8036     0.0254    0.0204   1958.5399    0.9991   ⋯
        w[8]     1.7987    0.8424     0.0266    0.0165   1915.4226    0.9994   ⋯
        w[9]    -0.0478    0.7916     0.0250    0.0188   2057.6024    0.9993   ⋯
       w[10]     1.9124    0.7558     0.0239    0.0177   1717.8319    1.0001   ⋯
       w[11]     0.6770    0.8081     0.0256    0.0194   1902.8984    0.9991   ⋯
       w[12]     1.2699    0.7636     0.0241    0.0167   1858.3391    0.9991   ⋯
       w[13]     1.2766    0.7838     0.0248    0.0143   2455.5405    0.9990   ⋯
       w[14]    -0.7512    0.7893     0.0250    0.0163   1702.2662    0.9990   ⋯
       w[15]    -2.8963    0.8600     0.0272    0.0146   1818.3626    0.9993   ⋯
       w[16]    -2.7241    0.9117     0.0288    0.0206   2450.0885    0.9992   ⋯
           a   104.2644    0.3407     0.0108    0.0131    952.1175    1.0038   ⋯
           σ     6.0974    0.1558     0.0049    0.0031   1775.3299    0.9990   ⋯
                                                                1 column omitted

Quantiles
  parameters       2.5%      25.0%      50.0%      75.0%      97.5% 
      Symbol    Float64    Float64    Float64    Float64    Float64 

        w[1]    -2.5281    -1.2683    -0.7015    -0.1222     1.1072
        w[2]    -2.5470    -1.3879    -0.8532    -0.2840     0.8402
        w[3]    -0.6997     0.2898     0.8241     1.3810     2.3930
        w[4]    -0.6385     0.2026     0.6853     1.1511     2.0438
        w[5]    -1.0201     0.0450     0.6056     1.1317     2.2009
        w[6]    -2.6198    -1.7510    -1.2338    -0.7059     0.1801
        w[7]    -1.3066    -0.3101     0.2989     0.8075     1.7683
        w[8]     0.1302     1.2344     1.8055     2.3861     3.4671
        w[9]    -1.6704    -0.5649    -0.0236     0.4951     1.4676
       w[10]     0.4834     1.3700     1.9061     2.4198     3.3756
       w[11]    -0.8996     0.1316     0.6925     1.2162     2.2533
       w[12]    -0.1711     0.7374     1.2467     1.7692     2.7470
       w[13]    -0.1589     0.7075     1.2470     1.8217     2.8017
       w[14]    -2.2540    -1.2961    -0.7754    -0.1916     0.7812
       w[15]    -4.5643    -3.4532    -2.9137    -2.2586    -1.3194
       w[16]    -4.3992    -3.2736    -2.7774    -2.1793    -0.8404
           a   103.5640   104.0427   104.2656   104.4941   104.9266
           σ     5.8152     5.9908     6.0953     6.1964     6.4151